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The dynamical evolution of self-gravitating scalar field configurations in numerical relativity is 
studied. The previous analysis on ground state boson stars of non-interacting fields is extended to 
excited states and to fields with self couplings. 

Self couplings can significantly change the physical dimensions of boson stars, making them much 
more astrophysically interesting (e.g., having mass of order O.f solar mass). The stable (5*) and un- 
stable (U) branches of equilibrium configurations of boson stars of self-interacting fields are studied; 
their behavior under perturbations and their quasi-normal oscillation frequencies are determined 
and compared to the non-interacting case. 

Excited states of boson stars with and without self-couplings are studied and compared. Excited 
states also have equilibrium configurations with S and U branch structures; both branches are 
intrinsically unstable under a generic perturbation but have very different instability time scales. 
We carried out a detailed study of the instability time scales of these configurations. It is found 
that highly excited states spontaneously decay through a cascade of intermediate states similar to 
atomic transitions. 

PACS number(s): 04.40.-b, 04.40. Dg 



I. INTRODUCTION 

Various particle physics models suggest that bosons might play an important role in the evolution of the universe. 
These models predict the abundant production of these bosonic particles in the early universe and their presence 
in large numbers even today. Although the bosonic particles have never actually been directly detected they are 
considered as leading candidates of dark matter [1]. These bosons could by a Jeans instability mechanism [2] condense 
into compact gravitating objects like boson stars. 

Boson stars are made up of self-gravitating complex scalar fields with or without further self coupling [3, 4]. The 
equilibrium configurations represent an exact balance between the attractive effect of gravity and the natural tendency 
for the scalar field to disperse. The stability of such an object is hence a central issue. It is well known that equilibrium 
configurations of boson stars have stable (S) and unstable (U) branches as well as a hierarchy of ground and excited 
states [3] -[13]. In the first paper in this series [5], the dynamical evolution under various perturbations of ground state 
boson stars made up of non self-interacting scalar fields was studied. 

In the absence of self coupling the mass profile of the ground state configurations, when plotted against the central 
density 0(0), has a peak at M = 0.633mp;/m (where mpi refers to the Planck mass) corresponding to a mass of about 
10^^ kg at (j){0)c = 0.271 (for bosons of mass m = IGeV). Configurations for which (j){0) < (/'(O)c (the S branch) are 
stable to perturbations [5,12] while those with (j){0) > (j>{0)c (the U branch) are unstable. Stability here refers to the 
ability of a 0(0) < 0(O)c star to settle to a new configuration in the same branch under perturbations. A 0(0) > (f>{0)c 
configuration star is unstable in that, upon perturbations, it cannot stay on the same branch. If it cannot lose enough 
mass and settle to a stable state it either collapses to a black hole or disperses to infinity. Stable boson stars have very 
specific quasinormal modes of oscillation under perturbations, a feature important for the detection and identification 
of these stars. 
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In the present paper the stndy of paper I[5] is extended to the excited state boson star, and to the case of a A^* 
self-interacting scalar field. The dynamical evolutions of such systems are studied numerically. 
The action for the system studied in this paper is given by 



/ = 



1 



(1.1) 



There are two reasons to include the self-coupling interaction[6,8]. Firstly without the A^* term, the maximum mass 
of a boson star, 



M = 0.633mp;/m 



(1.2) 



(where mpi refers to the Planck mass and m the mass of the boson) could be too small to be astrophysically significant. 
For example, for m =lGeV, M = 10~^^Mq where Mq is the solar mass. On the other hand, for interacting fields, 
even with a small coupling, the mass of the star can be large [6]. In this case 



M ■ 



0.06%/A^, 



(1.3) 



which is larger than (1.2) by a factor of ~ V^^^^ ~ 10^^-\/A, for m = IGeV. A moderate value of A = 0.01 then leads 
to M = O.IMq, which is particularly interesting due to the gravitational microlensing data[14]. 

Secondly boson stars give us a way to study local anisotropy and its effects. The larger the self coupling parameter A 
the smaller the fractional anisotropy. Changing the self coupling parameter for a given central density provides a way 
to vary this anisotropy in a natural way [8]. By anisotropy we mean that the radial and tangential components of the 
pressure are different. This is of interest becaiisc deviations from perfect fluid assumptions for even nuclear matter is 
expected in the presence of strong gravitational fields. In boson stars this anisotropy appears very naturally. Although 
we have not studied specifically the anisotropics in this paper, it provides motivation for adding a self-coupling term. 

The results of Paper I are also extended by considering the evolution of excited states with and without self coupling. 
This is of importance because if boson stars exist and are detected, they are most likely those interacting with their 
environment and going through some excitation process. Excited configurations might also be intermediate stages 
during the formation process of these stars. In this study the configurations considered are spherically symmetric; all 
perturbations of equilibrium configurations are purely radial. Full 3D simulations are underway and will be reported 
in a future paper. The remainder of this paper is organized as follows: 

Section II sets up the mathematical foundations of the problem, including the equilibrium and evolution equations. 
The calculation of the initial data sets are discussed and the techniques used to evolve the system numerically are 
briefly outlined. 

Section III details the evolution of ground state configurations with self coupling. We show that they have similar S 
and U branch structures as boson stars without self-coupling. S branch stars are stable with regard to perturbations. 
By this we mean that under small perturbations they return to configurations on the same branch although not to the 
same configurations. We studied in detail the fundamental quasinormal modes of oscillation of S'-branch stars (those 
for which (f>{0) < (f>{0)c)- They are important characteristics for observations. They can also be used to predict the end 
point of evolution of perturbed stars, as well as a comparison between modes for different self couplings. In the next 
subsection of section III, the migration of a [/-branch star (that for which (j){0) > (j){0)c) to the 5-branch is described. 
The essential features of the A = case are retained. If a [/-branch star is perturbed by the addition of mass, the 
star will collapse to a black hole. When, as a result of perturbation, the star's mass is reduced, corresponding to 
annihilation of scalar particles, the star expands and moves to the 5-branch, oscillates, and settles to a new equilibrium 
configuration of lower mass. 

Section IV studies various aspects of the evolution of excited states. Excited states have similar band structures to 
ground state stars. In this paper, we study generic perturbations that may exist for a boson star in an astrophysical 
environment, e.g., some additional scalar particles falling in. The S-branch excited states have previously been found 
to be stable under infinitesimal perturbations that conserve the total mass and particle number of the boson star [4]. 
We find that these stars are inherently unstable irrespective of whether they lie on the S or the U branch but the 
time scales of instability are different. This result is consistent with the study of infinitesimal perturbations [12]. If 
they cannot lose enough mass to transit to the ground state they either collapse to black holes or, as in the case of 
stars for which M > N m (M = mass of the star, N = number of bosons and m = mass of one boson) disperse to 
infinity. The decay of some higher excited configurations are also studied. These higher node configurations cascade 
through intermediate configurations of lower excited states on their way to collapse. This is reminiscent of atomic 
transitions where atoms go from an excited state to lower states through intermediate ones, lending credence to the 
idea that boson stars are like gravitational atoms [3]. A brief conclusion follows in Section V. An appendix at the end 
shows some features of the the high A (A = X/Airm^G) configurations, including a calculation of quasinormal modes. 
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II. FORMULATION AND EQUILIBRIUM MODELS 



In this section the mathematical formulation of the problem, the creation of the equilibrium and perturbed boson 

star modes and the numerical code used to study them are described. The formulation is the same as that of Paper I in 
this series [5] except for the self- interaction term. The numerical treatment used in this code has various improvements 
over that described in Paper I. Some details of the numerical code, e.g., convergence tests etc, which are similar to 
those reported in Paper I will not be repeated here. 

The action for a self-gravitating scalar field given by equation (1), leads to the scalar field equation 



0, 

for the complex scalar field $ = $i -|- and the Einstein field equations 



(2.1) 



R 



SttGT, 



The metric for this spherically symmetric system can be written as 



(2.2) 



where g, the radial metric, and N the lapse are functions of (t,r) with r being the circumferential radius. This 
form of the metric is known as the radial gauge. In the absence of a shift vector this form of the metric can 
be maintained for all time by enforcing the polar slicing condition. This is a condition on the lapse that requires 
Keo + K^^ = where Kij is the extrinsic curvature tensor. This slicing condition causes the lapse N to decrease 
rapidly if an apparent horizon is approached [15]. 

The equilibrium boson star configurations are those in which the metric is time independent. The scalar field $ 
itself oscillates with fixed frequency wq: 

$(t,r) = $o(r)e-^-»t, (2.3) 

but due to the U{1) symmetry of the Lagrangian, the stress energy tensor and the spacetime geometry are time 
independent. In dimensionless coordinates we have 

r = mr, < = cjot, a = \/47rG$, A/' = N— , A: 



and the Einstein and Klein-Gordon equations under these conditions are 

o-Q = Xi 



m247rG' 



(2.4) 



(2.5) 



X'l 



1 9' 
- + — 
r r 



2 2 



Xi 



1 

iV2 



1 



(2.6) 



9 5" 2 3 

+(^org 

r r 



1 

iV2 



(2.7) 



N' = - 
2 



N Ng"^ 



2 2 

N 



(1 _ iv2) + rNxl - 2 Ag'iVra^ 



(2.8) 



where co = <1>oV47tG. A prime denotes d/dr and an overdot denotes d/dt. All quantities in this paper are reported 
in terms of these dimensionless parameters unless explicitly stated otherwise. Regularity at the origin requires that 
g{r = 0) = 1 and that all other quantities are finite at r = 0. For the solution to represent an isolated star, it 
is required that a{r = oo) = 0. This constitutes an eigenvalue problem. For each choice of a{r — 0), the above 
set of equations has a solution only when N{r = 0) takes on certain values. Different eigenvalues correspond to a 
different number of nodes in the solution of cr(r). Solutions are also obtained for different values of the coupling 
parameter A. Different families of equilibrium configurations are shown in Fig. 1. The mass profile of ground state 
boson stars with self coupling has the familiar structure seen in non self- interacting fields [5] which is also found in 
white dwarfs and neutron stars. The mass grows to a maximum as the central density is increased, and then decreases 
with further increase in central density. (See Fig. 1). The maximum mass increases with A but the profile is similar. 
The expectation that the branch to the left of the maximum is stable (5-branch) while that to the right (?7-branch) is 
unstable, as in the case of groimd state configurations without self coupling, is found to be true. Stability here refers 
to the ability of S'-branch stars to settle to new iS'-branch configurations when perturbed, a feature of importance for 
the long term existence of these stars. 
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A. Evolution equations 



The configurations described above are time independent, equilibrium solutions to the Einstein equations. The aim 

of this paper is to study their dynamical properties according to the coupled Einstein- Klein-Gordon equations. In 
this paper only spherically symmetric configurations are studied. In actual numerical evolution the following set of 
variables are chosen 



ipi = rai , ip2 



where 



1 dTpi 

rC72, TTl = 

a at 



N 

a=—, 
9 



_ 1 9^2 
7r2 = KT^ 

a at 



(2.9) 



(2.10) 



and the subscripts on ipi denote the real and imaginary parts of the scalar field multiplied by r. 

In terms of these variables and the dimcnsionless ones in the previous section the evolution equations are as follows: 
The radial metric function g evolves according to 



g = N{TTia[ + 7r2(T2)- 

The polar slicing equation, which is integrated on each time slice, is given by 



2 



-+r[(CTi) +(ct2) -5(^1+0-2)]+ ^(«^l+'72) 



The Klein-Gordon equation for the scalar field can be written as 



(2.11) 



(2.12) 



1,2, 



(2.13) 



Ipi = awi, i = l,2. 
The hamiltonian constraint equation is given by 



2g' 0^-1 tt^+ttI a[^+a'i 



A. 



2\2 



J-2g2 



K + ^2) -2(^1+^2) 



0. 



(2.14) 



(2.15) 



It is not solved during the evolution but as it is in principle conserved by the evolution equations, it is monitored 
closely as an indicator of the numerical accuracy of the simulation. For further details see [5]. 



B. Boundeiry conditions 

Regularity conditions require that g{r = 0) = 1, and g,N,ai and (T2 have vanishing first spatial derivatives at 
r = 0. To implement this condition numerically, the range of r is extended to include negative values; g, N, ai and 
(72 are required to be symmetric about r = 0. In addition tpi, tp2, tti and 772 are antisymmetric about r = 0. The 
antisymmetry allows the determination of ^, at the origin as the first derivatives of ipi at r = 0. The value of the lapse 
function is fixed at the outer edge on each time slice. Its value at the origin is determined by integrating Eq. (2.15) 
inward from the outer boundary. The value of g is determined by the evolution. The mass of the star is determined 
by the value of g at the edge of the grid: 



,^ 1 
M = -r 
2 



1 



5f2(oo) 



'Pi 



m 



(2.16) 



where mpi is the Planck mass and m is the mass of the boson making up the star. 

The boundary condition on the scalar field is an outgoing scalar wave condition. However, since the dispersion 
relation of the massive scalar field is non trivial, 

a^k^ = uj^ - N'^rr?, (2.17) 

(where a = N / g) there is no perfect algorithm for the implementation of the outgoing wave condition. Here we have 
adopted a two tier approach: 
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• i) A "sponge" regioii[5] is constructed by adding a potential term at the outer edge of the computational domain. 

tn — D < r < rN 

where rpf is the r value of the outermost grid point and D is an adjustable parameter representing the width 
of the sponge. D is typically chosen to be a few times the wavelength of the scalar radiation moving out. The 
extra potential term in the above equation is designed to allow waves to propagate outward but damp incoming 
waves. 

•ii) 

At the outermost grid point we require 

i/j = -aip' - —tp. (2.19) 

This is an exact outgoing wave condition only in the case m = 0. The second term on the right hand side 
represents the finite m correction to leading order (recall N = N^). The sponge is designed to absorb the 
reflection coming from this approximate outgoing wave condition, note that recent work on this problem 
following a hyperbolic approach seems to provide a simple and more accurate outgoing boundary condition [16]. 

III. DYNAMICAL EVOLUTIONS OF PERTURBED GROUND STATE STARS WITH 

SELF-INTERACTION 

A. Nature of perturbations 

We study the dynamical properties of the boson stars by perturbing the equilibrium field distribution. The accretion 
or annihilation of scalar particles is simulated by the addition of field in the outer regions of the star or by decreasing 
it in denser regions of the star respectively. Another type of perturbation that has been effected is changing ipi and 1^2 
of the equilibrium configuration. This perturbation changes the kinetic energy density distribution. In either case the 
changes in the metric functions g and N are determined by the constraint equations and the polar-slicing condition 
(integrating Eqs. (2.12) and (2.15) on the initial slice.) The magnitude and the length scale of the perturbations can 
be chosen arbitrarily. The perturbations are always spherically symmetric. 

In the rest of this section we present results obtained on the dynamical evolutions of ground state equilibrium 
configurations perturbed in this manner. Evolutions of both ^-branch and f -branch stars with and without self 
coupling are considered. 

B. S-branch perturbations 

As shown in Ref. [5], in the free field (A = 0) case, a perturbed 5- branch star oscillates with a definite frequency, 
losing mass through bursts of scalar radiation at each expansion and finally settles to a new S'-branch configuration 
of lower mass. Here the effect of a self-interacting A term on this behavior is examined. In Fig. 2 a typical example of 
the perturbed field configuration and radial metric of a star of A = 10 with a central density of cr(0) =0.1 is shown. 
This star has been perturbed by accretion of scalar particles in a region of lower density. Its evolution is detailed 
below. 

Fig. 3a shows the radial metric as a function of distance for the same configuration at various times. The labels 
A, B, C, D correspond to times (in units of the inverse of the underlying scalar field frequency) t = 192, 306, 391, 505 
respectively. The positions of the peaks are labeled Rq,Ra,Rb,Rc,Rd, where Rq is the position of the initial 
unperturbed peak. Here Rq = 7.95, Ra = 8.55, Rb = 6.6, Rc = 8.2 and Rn = 6.65 where the length scales are in 
terms of the inverse mass of the boson. The oscillations are shown clearly in Fig. 3b which is a plot of the maximum 
value of the radial metric as a function of time. The point where this function is a maximum corresponds to the 
core of the star contracting to its minimum size in a cycle. Similarly, the maximum radial metric starts to decrease 
as the star expands. At each expansion the star loses mass through scalar radiation. The oscillations damp out in 
time as the star starts settling to the new configuration. A plot of mass vs. time is shown in Fig. 3c. The amount of 



gN- 



1,2, 



(2.18) 
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scalar radiation decreases in time as the oscillations damp out, as can be seen from the figure. The slope of the curve 
steadily decreases as the star starts settling down to its new lower mass configuration. The mass is measured at the 
inner edge of the sponge. Exact details of the curves have some dependence on the sponge parameters but the basic 
results are the same. 

A characteristic of the boson star that could be important for its observation and identification is its fundamental 
oscillation frequency which can be determined from Fig. 3b. We found / = 1/(199 A'' (oo)) =4.7x10^'^. The oscillation 
frequencies for a large number of S-branch stars have been compiled in this way. Fig. 4a shows a plot of the oscillation 
frequency versus mass for many slightly perturbed configurations (masses within 0.1% of the unperturbed mass.) As 
the mass increases, the frequency increases and then drops down as the transition point ( fJ^Q^ = 0) is approached, 
signaling the onset of instability. This is seen for both the non-self interacting as well as the self-interacting case. 
These quasinormal modes of oscillation characterize 5-branch stars. The point of transition from the S to the U 
branch corresponds to a zero frequency of oscillation [12]. 

For a given mass higher A S'-branch stars have a lower oscillation frequency than similar mass lower A stars, unless 
one is near the transition point of the lower A configuration. This is not too surprising, since for a given mass the 
radius of the star increases with increasing A. We have seen this trend even for A values as high as 1600 (see discussion 
on high A stars in the Appendix) . However since the maximum mass of higher A configurations is greater than that 
of lower A configurations their maximum oscillation frequency could be greater than that for lower A stars. (This 
can be understood as a size effect: on the S-branch, higher mass stars are smaller and have higher frequencies.) This 
can be seen in Fig. 4a. The maximum oscillation frequencies of A = 5 and A = 10 configurations are higher than 
that of the A = case. As the stars get much larger, though, the highest frequency starts to decrease. For example 
the maximum frequency for A = 30 stars being less than that of A = 15 stars which is in itself lower than that of 
the A = 10 case. A perturbation calculation for the high A case is shown in the appendix to show the dependence 
of quasinormal mode frequencies on A for high A configurations. The frequency is proportional to the inverse of the 
square root of A. Thus as the stars get really big they oscillate less and less rapidly and numerically it is no longer 
feasible to evolve them. (The time step used in the numerical simulation cannot be increased as it is determined by 
the intrinsic oscillation time scale of the scalar field, which is many orders of magniture shorter than the oscillation 
time scale of the whole star for these cases.) 

The quasinormal mode curves are also useful in determining evolutions of strongly perturbed S-branch stars and 
the final configurations they could settle into. A perturbed star loses mass and settles to a final configuration 
corresponding to a position on the solid line in the figure. In Fig. 4b we single out the A = 10 curve and plot the 
evolution of the S-branch star discussed above. The points PI, P2, P3 show the route to a new configuration. These 
points correspond to times t = 0, 1200 and 4800 respectively. By extrapolating this line to where it meets the A = 10 
curve one could expect a final mass of about .76 mp;/m . 

C. ?7-branch perturbations 

For the case of non-self interacting fields it was shown in Ref. [5] that accretion of scalar fields causes tZ-branch 
stars to collapse to black holes. However, lowering the density can make the star migrate to the S-branch. These 
features are also seen in the presence of self coupling. Fig. 5 shows a migrating A = 30 star whose unperturbed overall 
field density a has been decreased by about 10%. Fig. 5a shows the behavior of the radial metric in time as a function 
of radius. It oscillates about the final S-branch configuration that it will settle into. This final state is shown on the 
plot as a dark line. Fig. 5b shows the maximum radial metric as a function of time. The star initially expands rapidly 
as it moves to the S-branch. This can be seen from the sharp drop in the radial metric. Once it moves to the stable 
branch, it oscillates about the new configuration that it is going to settle to. Fig. 5c shows the mass of the star as 
a function of time. It loses mass at each expansion losing less and less mass at each subsequent expansion and the 
curve gets smoother and smoother as it prepares to settle to its final state. 

Fig. 6, which shows the oscillation frequency as a function of mass for A = 30, can be used to predict the end point 
of migration. Points Ql, Q2, Q2> and (54 show the migration of this star. These correspond to times of 500, 1000, 2000 
and 3500. The oscillation is clearly damping out. The final configuration it is expected to settle down to is shown as 
a dot and corresponds to a stable star of central density (t(0) = 0.0817 with a mass of 1.037mp;/m. This example is 
typical of a number of simulations of [/-branch ground state configurations with self coupling. 

For higher central density stars on the JJ-branch, the mass versus central density curve has a second, gentler peak, 
similar to the white dwarf neutron star situation. One might suspect that this corresponds to another stable and 
unstable branch respectively. However we find that configurations on both sides of the peak are unstable. These 
configurations always disperse upon perturbation, consistent with the fact that they have M > N m, where M was 
the mass of the star and A'' m was the number of bosons multiplied by the mass of a boson. 
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In summary ground state configurations of boson stars with self coupling have stable {S) and unstable (U) branches 
just like boson stars without self coupling. The stable configurations have very specific quasinormal modes of oscil- 
lation. The addition of a self-coupling term serves to increase their mass. The cases considered so far correspond to 
A ~ 10^ or less. Appendix A deals with very high A stars which need a different technique due to the very different 
timescales involved in these evolutions. 



IV. EVOLUTIONS OF EXCITED STATES 

Excited states of boson stars have field configurations characterized by nodes. The first excited state has one node, 
the second has two and so on. Studying their stability and the time scale of decay is important in determining the 
likelihood of finding them in nature. 

The mass profiles of excited state boson stars are similar to ground state stars. Fig. 7 shows the mass versus central 
density curves for ground, first and second excited states of boson stars without self coupling. The maximum mass 
increases with the number of nodes as expected. The similarity of mass profiles of excited boson stars to their ground 
state counterparts, might lead one to expect stable and unstable configurations to the left and right of the maximum 
mass respectively in analogy with ground state configurations. However, our numerical studies show that the excited 
boson star configurations on both sides of the peak are inherently unstable except that the time scales for instability 
are different. If they cannot lose enough mass to go to the ground state they become black holes or totally disperse. 
This occurs even if no explicit perturbations are put in the numerical evolution other than those introduced by the 
finite differencing error in the numerical integration. We have also carried out perturbations that correspond to more 
scalar particles falling on to the star, or those that decrease the scalar field strength at the center point corresponding 
to scalar particles decaying through some channels [5]. The instability shows up in all cases studied. We note that 
this instability is NOT in contradiction with the result of [4] , which concluded that S branch excited states are stable 
under infinitesimal perturbations that strictly conserve M and A'', (where M is the mass of the star, A'' the number 
of bosons and m the boson mass.) Our result of instability under generic perturbation is consistent with the studies 
of A = stars under infinitesimal perturbation[12]. The presence of a self coupling term increases the time scale of 
instability but the essential pattern remains the same. In the following we will first give a detailed account of the 
dynamical evolution and carry out a study of the instability time scales. 

Excited state stars with masses less than the maximum mass of a ground state configuration always form ground 
state stable configurations. In fact even stars with masses somewhat greater than this can lose enough mass during 
their evolution so as to go to the ground state. Fig. 8a shows a perturbed 1-node star whose mass has been reduced 
by about 8% to 0.9 rripi/m by a perturbation making a transition to the ground state although the mass is greater 
than the maximum groimd state mass of 0.633mp;/m. A substantial amount of scalar radiation is emitted in the 
dynamical evolution, which brings the mass below the critical value. The evolution of the radial metric function is 
shown. Although the plot is shown only to a radius r = 100, the actual evolution was carried out to r = 300. The 
two peaks at t = is indicative of a first excited state. One of the peaks disappears gradually as the star goes to the 
ground state. The star then oscillates about the ground state configuration that it will finally settle into. In Fig. 8b 
we show a 3 node configuration with a total mass of 0.92top,/to, going to the ground state after radiation by scalar 
waves carries off the excess mass and kinetic energy. We have plotted the density function against the radius and 
time of evolution. By the density function we mean density p multiplied by an factor which is the mass per dr at 
radius r. pr^ has n + 1 maxima for an n node star and hence here we have 4 sets of lines initially. At the end of the 
simulation, we see that it settles down to a ground state configuration with small oscillations with ever decreasing 
amplitude. For these simulations of low central density stars we put in an explicit perturbation to the equilibrium 
configuration since the instability time scales are extremely long without that. 

For stars with higher central density, there is a critical density above which the stars cannot lose enough mass to go 
to the ground state but collapse to black holes. In our numerical simulation for one node stars this critical density is 
(t(0) = (72 = 0.048. As the central density is increased the kinetic energy of the highly compressed initial equilibrium 
configuration is increased and the star first expands before the eventual collapse to a black hole. As the central density 
further increases towards the M = Nm point the expansion phase becomes longer. In Fig. 9a we show the density 
fimction (p x r^) against radius at various times for 4 configurations (an S branch configuration with cr(0) = 0.1 and 
three U branch configurations of central densities (t(0) = 0.3, cr(0) = 0.4 and cr(0) = 0.5). The initial configurations 
are the equilibrium ones without any explicit perturbation (except those introduced by the discretization used in 
the numerical simulations). The first frame shows the S branch star radiating a little as it makes a transition to 
the ground state. However it cannot sustain this state for long and it collapses to a black hole. The time of decay 
decreases in the case of a [/ branch stars of central density 0.3 and 0.4. However for the cr(0) = 0.5 star the star 
is clearly more dispersive than the previous ones. It goes through an expansion phase initially although it finally 
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collapses to a black hole. Configurations with a{0) > 0.54 have M > Nm. They do not collapse to black holes but 
disperse to infinity. 

Next we turn to a study of instability time scales. In simulations where a black hole will form, the imminent 
development of an apparent horizon leads to a rapid collapse of the lapse due to the polar slicing used in the evolution. 
We take the time of collapse of the lapse at the origin to ~ 10~^ of its initial central value to be the approximate 
time for formation of the black hole. Fig. 9b shows this time scale for a 1-nodc star without self coupling. We plot 
the decay time scales of first excited state stars as a function of central field density. Again no explicit perturbation is 
applied in these evolutions other than the discretization error in the simulation. In order to make a fair comparison 
of the time scale due to such a perturbation we cover the radius of the star in all cases by the same number of grid 
points. The maximum ground state mass for stars without self coupling is around .633 mp^/m. This corresponds to a 
central density of cti = 0.021 for a 1-node star. We described earlier that stars with central densities below a2 = 0.048 
(mass of 0.91mp;/m) lose enough mass and move to the ground state. Beyond that and upto the M = N m, point they 
collapse to black holes, while stars with M > Nm corresponding to a central field density of (t(0) > 0.541 disperse to 
infinity. The time by which collapse takes place to a black hole decreases with increasing central density along the S 
branch. This trend continues for a while into the U branch (starting at cr = 0.25 until a ~ 0.4) but as one approaches 
the M = Nm point the stars lose a significant amount of matter to infinity before they collapse to black holes and 
evolve on a longer time scale. For example a star of central density 0.5 has an initial radius of r ^ 9. (The radius of 
the star is defined as the radius which contains 95% of the mass of the star.) Its radius increases to as much as 115, 
more than an order of magnitude, before it starts collapsing. Dispersion time scales of a couple of stars for M > Nm 
which disperse to infinity (instead of collapsing to a black hole) are also shown on the figure. To give a sense of the 
instability time scales of these configurations we take the time scale to be the time by which these stars disperse to 
ten times their original radii. This time drops drastically for stars with M » Nm. 

Next we turn to the case of A 7^ 0. Fig. 10a shows black hole formation for a A = 30 star in the first excited state, 
with a central density (t(0) = 0.1. This star had an initial radius of about 20.7 where the radius is again defined as 
that containing 95% of the mass. The lapse finally collapses to zero indicating that an apparent horizon is about to 
form. The time scale of collapse to the black hole was around 1985 compared to a time scale of less than 800 for a 
cr(0) = 0.122 star of similar radius without self coupling which is shown in Fig. 10b. This is to be expected as the A 
term represents a repulsive force. The time is again determined by the lapse collapsing to 10~^ of its original value 
at r = 0. 

We now turn to the evolution of a highly excited state. In Fig. 11a we show the initial field configuration of a 
star containing 5 nodes. For a five node star the density has a central maximum, and then five local maxima, each 
subsequent one smaller than the one preceding it. This star is then evolved without any perturbation except those 
introduced by the discretization error of the numerical evolution. In Fig. lib we show a contour plot of the evolution 
of the density function in time, pr"^ has n + 1 maxima for an n node star and hence here we initially have 6 sets of 
lines centered at dimensionless r = 5, 13, 35, 52 and 75 respectively. This star has central density cr(0) = 0.075 and 
it collapses to a black hole after a long evolution. In the process we see intermediate states with fewer numbers of 
nodes. For comparison in Fig. 11c show the contour plot for the 5 node star upto a time of t = 1000 at which time 
it has decayed into a four node state, against the equilibrium density function of a 4 node star with central density 
a = 0.06. Very clearly the max;ima are at similar radii (although the size of the peaks are somewhat difli'crcnt). 

This feature of nodes disappearing and the star cascading through lower excited states is characteristic of decay of 
higher excited states of boson stars. Although in this case and at this time the decaying star is close to a specific 
lower excited state, in general the decaying star is roughly a combination of lower excited configurations. This is 
similar to the decay of atoms in excited states. However, we note that for the "gravitational atom" [3] there is no 
exact superposition due to the intrinsic non-linearity of the system. 

V. CONCLUSION 

In the first paper in this series, the behavior of boson star ground state configurations under various perturbations 
were reported. In this paper the study has been extended to include boson stars of self-interacting field and also the 
behavior of boson stars in the excited states. 

The self-coupling term is important as it can have dramatic effects on the the mass of the boson stars [6] leading to 
boson stars of the order of a solar mass. The mass profile retains the features of boson stars without self-coupling, 
having a central maximum with a stable branch and an unstable branch. All configurations to the left of the central 
maximum in the mass vs. central density curve (see Fig. 1) are stable. Under small perturbations they have very 
specific quasinormal modes of oscillation as seen in Figs. 4 and 6 and under perturbations they settle down to new 
configurations on the same branch. Configurations that lie on the unstable or i7-branch, either migrate to new 
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configurations on the 5'-branch or collapse to black holes, when perturbed. These are characteristics shared by boson 
stars with or without self-coupling . 

Excited states arc configurations with nodes. The field of an n*'' excited state star has n nodes and its radial metric 
has n+ 1 peaks. Their mass profiles are similar to the profiles of boson stars in the ground state, which makes it appear 
as if they have a stable and an unstable branch of configurations. However, irrespective of which branch they lie on, 
excited boson stars are unstable with different instability time scales. Low density excited stars having masses close 
to ground state configurations will form ground state boson stars after evolution. Denser configurations form black 
holes with the decay time decreasing with increasing central density till one approaches the density corresponding 
to zero binding energy. As the central density approaches this central density the kinetic energy of the star starts 
to increase as it becomes more dispersive. It still collapses to a black hole but on a larger time scale. Beyond this 
point for densities corresponding to positive binding energy the stars disperse to infinity. We studied the time scales 
of their instability in figure 9. 

An interesting feature in the collapse of excited state boson stars is that they cascade through intermediate states, 
during this process, rather like atoms transiting from excited states to the ground state, suggesting that boson stars 
behave in some ways like gravitational atoms [3]. However, an investigation of the possible decay channels (selection 
rules) seems much more difficult (if at all possible or meaningful) here, due to the intrinsic nonlinearity of the theory. In 
this paper we have reported evolutions of spherically symmetric configurations. We are currently extending the study 
to full 3-D without spherical symmetry. The numerical study of 3D boson stars in addition to being an interesting 
physical problem is also a testbed for 3-dimensional numerical codes, that enable us to study compact self-gravitating 
objects without having to deal with hydrodynamic sources as in neutron stars, and singularities as in the case of 
black holes. In particular we aim to study the general two body problem in relativity by evolving two 3-dimensional 
scalar field configurations. The inspiral coalescence of such systems could have interesting physical implications as 
the gravitational wave emitted does not sensitively depend on the internal structure of the compact objects until the 
late stages of coalescence. Studying the 1-d behavior has been an important tool in testing our 3-d codes, providing 
testbeds in the spherically symmetric limit. 
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APPENDIX A: HIGH A CASE- 



EQUILIBRIUM CONFIGURATION AND QUASINORMAL MODE 
DETERMINATION 



While calculating the eigenvalue for the equilibrium boson star it is found that it gets increasingly difficult to 
calculate the eigenvalue as the value of A gets large, because the equations get very stiff. There are two scales to the 
problem: a scale of slow variation of the field inside a certain radius related to A followed by rapid decay outside it. 
This makes an effective surface layer to the star making it more similar to neutron stars. It turns out that the large 
A limit can be treated using a set of approximate equations that are exact in the A = oo limit [6]. By making the 
following change of variables: 



r = r/VA., G = \fh.a 



(Al) 



the equilibrium equations reduce to 
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whore the primes refer to differentiation with respect to r. In the Hmit of A = oo, one can keep terms to leading order 
in in (A2)-(A4). In particular (A2) reduces to 

N={a^ + l)-i. (A5) 

However, we note that this is valid only for a ground state configuration. For a state with nodes it would not be 
reasonable to neglect derivative terms compared to terms proportional to a, which is zero at a node. 

To get an estimate of the accuracy of the high A approximation, we compare the solution using the approximate 
equation for high A to the brute force numerical solution of the complete set of equations, for a A = 800, a = 0.05 
boson star in the ground state in Fig. 12. The agreement of the fields is quite good until the outer region where 
the approximate equations cause the field to abruptly fall to zero. Comparisons of the radial metric and the lapse 
arc also shown. In the A = oo limit the approximate equations arc exact and the star really has an outer surface 
reminiscent of a neutron star. In fact, an equation of state can be written[6]. In Fig. 13 the mass and particle number 
versus central density a for high A stars is shown. One expects as in the case of other ground state configurations 
that the configurations with M > NpUi disperse when perturbed while those on the U branch with M < Npun would 
be unstable and if unable to migrate to the S branch under perturbations would form black holes. (Here we use the 
symbol Np for the particle number and not N so as not to confuse it with the lapse as both these functions figure 
prominently in the analysis that follows.) The particle number is calculated from the current J° and is given by 

Np^ A-K j r'^j^ua^dr. (A6) 

Here a has been replaced by a; cr and d?r by Anr^dr for the spherically symmetric case. In terms of the bar coordinates 

we see that Np — ^/ANp. 

Next we turn to the determination of the quasi-normal frequency (QNM) of the high A stars. In principle one 
could determine the QNM using the dynamical studies as performed for the A = case. However the procedure is 
extremely computationally expensive. The scalar field has an inherent oscillation of about 27r and the evolution time 
steps must be small enough to resolve it. However the code must run long enough to see a few metric oscillations 
in order to determine the quasinormal mode. As A gets large, the sizes of the stars also get large leading to a lower 
frequency of oscillation. In order to determine the QNM we use instead the following perturbation analysis based 
on[8,ll] but using our notation for lapse, fields, time, radius and self-coupling as defined in section II. 

We write the perturbed fields as: 

CT=((Ti+W72)e''^*, g = go + 6g, N = No+5N (A7) 

where 

0-1 = ao{r){l+5ai{r,t)), (T2 = ao{r)6a2{r,t). (AS) 
the Klein-Gordon equation can be written as 



^(#-f)(^i+^2)-2^< 
The field cto satisfies the equilibrium equation. 



Expanding to first order perturbations using (A9) and (AlO) we get 

\r N g aoj (tq \ go^ A/o 



2go(l + ^Aao^ J dg - go'^Aao'^Sai = 0. 



10 



Prom 



and 



(where T^j, = 47rGT^y) we get we get the equations 



5"Va'^,/3-(l + iA||c7r)|H| 



6N' No'SN 



go^ No J rgo^No rgo'^N^' 



and 



6%' 



50 \\NoJ No'' Nogo r\No go J \ No'^go 



Sg 



1 f5N" No',^, , No 



+ 2 



50^ V ^0 

No 



No' 



5N' + 



^SN-^(s,'-^Sn 
iVo ^050 V .90 J Nogo 



fja' 



No 



No 



SN' - ^SN 
No 



1 f ,57V' No' 



r [ No No- 



go 50 



and the equations 



ST^ = -^6N - j^a'o'Sg + 5a, (^ao' + ^a'o' + 2(1 + iAao^)ao^) 



and 
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respectively. Adding 5Tq to 5Tl we get 
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Substituting equation (A20) as well as for 6a2 from (A16) and (A19) in (All) we get 
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Adding STq to 6T2 and substituting in equation (A20) and its derivative we get 
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Using the expression for the particle number Np = d^xJ^^Jg where J° = ig^^ {(t'fl^* — c.c.) we get 
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In terms of bar coordinates defined in the beginning of this section equation (A21) becomes 
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Similarly equation (A23) becomes 
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Using 5Np = (charge conservation) which is appropriate for large A where SNp is given by (A26) we get on putting 
this in (A24) 
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This equation suggests that Sg goes Ukc Writing Sg as —x^5g and 5a as —x^Sa wc then sec froni(A25) and 

(A27) that the quasinormal mode frequency x for high A configurations must go hke for a given a. We have 
numerically evolved stars with the same a and compared the QNM frequencies obtained to the inverse ratios of the 
square root of their A values for A = 800, A = 1200 and A = 1600 confirming the above analysis. In Table I we show 
the comparison between the perturbation analysis and the numerical result. We see that the perturbation result gets 
more accurate for increasing A. Wc note that configurations that have the size of a neutron star would have to have 
A of order 10^^ with QNM frequency of order x = 10~^^. 
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FIG. 1. The mass profiles of ground state boson stars for difi'erent values of self-coupling constant A are shown. The increase 
in mass with A is clear although the profiles are very similar. 
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FIG. 2. The comparison of a strongly perturbed ground state A = 10, S-branch star [M = 0.781m%i/m,a{0) = 0.1] to the 
unperturbed configuration [M = 0.722m%i/m] is shown. The solid lines correspond to the unperturbed configuration and the 
dashed ones to the perturbed star. The perturbation shown corresponds to the addition of scalar-field cr at i = 0. 

FIG. 3. (a) The evolution of the radial metric Qrr ~ for the configuration shown in Fig. 2. The initial perturbed 
configuration is labeled t = Q. The unperturbed configuration is also shown. The spatial distributions of the radial metric 
labeled A, B, C, D correspond to times t = 192, 306, 391, 505 respectively. The radial positions of the pealcs of the radial metric 
for these times are labeled Ra, Rb,Rc, and Rd- The values Rb and Rd are so close that they appear as one thick line in the 
figure, (b) The peak value of the perturbed radial metric is plotted over a long time. The points labelled A, B, C, D correspond 
to the same labels in (a). The oscillations decay in time, (c) The total mass of the star is plotted as a function of time. The 
mass loss through scalar radiation decreases in time as the oscillations start damping out. 

FIG. 4. (a)The oscillation frequencies of different ground state boson star configurations are plotted as functions of mass, for 
A = 0, 5, 10, 15, 30, 100 and 200. The curves are obtained by slightly perturbing (perturbed mass within 0.1% of the unperturbed 
mass) S'-branch stars. They reach a peak and then drop down at the approach of the majcimum mass allowed for a given A 
signalling a transition from stability to instability. The frequencies for a given mass for higher A stars are lower than those for 
lower A as a result of their increased size. However their overall maximum frequencies could get bigger than for lower A stars 
because of their increase in maximum mass. As can be seen A = 5 and A = 10 stars have higher maximum frequencies than 
do A = stars. As the stars get very much larger the maximum comes down as shown in the figure. A = 30 stars have lower 
maximum frequency than A = 15 stars which have lower maximal frequency of oscillation than do A = 10 stars. (b)The highly 
perturbed A = 10 S'-branch star of figure 3, has an oscillation frequency below the A = 10 solid line. Its movement towards the 
solid line is shown through points PI, P2 and P3 corresponding to times 0, 1200, and 4800 respectively. 

FIG. 5. (a) The radial metric = grr of a perturbed {/-branch ground state star is shown at various times. The curves 
1, 2, 3, a, b, c, d correspond to times t = 10, 20, 30, 340, 440, 540 and 640. The unperturbed star has a central density (t(0) = 0.23 
and a self coupling parameter A = 30. The initial equilibrium metric configuration has also been shown. The overall field 
density of this star has been lowered by about 10%. The t = curve corresponds to the initial perturbed radial metric. In 
the asymptotic region g^ is not oscillating but monotonically decreasing due to the mass loss, (b) The maximum value of the 
radial metric is plotted as a function of time. The initial sharp drop in the radial metric signifies the expansion of the star 
as it proceeds to the stable branch. There it oscillates about the new stable configuration that it is going to settle to. This 
corresponds to a star of mass M = 1.037mp;/m (whose metric configuration has been shown in figure 5a as a dark line). The 
points a and c correspond to two minima in the peak of grr which occur when the core of the star reaches its local maximum 
size. Likewise the maxima in the peak of grr at b and d correspond to the core of the star reaching its local minimum size, (c) 
The mass is plotted against time. The mass loss through scalar radiation at each expansion of the core (corresponding to the 
maximum radial metric reaching a minimum) decreases in time as the oscillations damp out. 

FIG. 6. The migration of the {/-branch star considered in the previous figure is shown after the star has moved to the 
S-branch. Points Ql, Q2, Q3 and correspond to times 500,1000, 2000 and 3500 respectively. Here too the mass loss 
decreases in time and the star finally settles to a stable configuration. 

FIG. 7. The masses of 0-node, 1-node and 2-node boson stars without self-coupling are plotted as a function of central 
density. The maximum mass of 1-node stars is 1.356mp;/m while the maximum mass for 2- node stars is expectedly greater at 
2.095mpj/m. The profiles are deceptively similar to their ground state counterpart. Excited state stars are inherently unstable 
irrespective of the branch they lie on, unlike ground state stars that can be termed stable or unstable depending on whether 
they lie on the branch to the left of the maximum mass or to the right respectively. 

FIG. 8. (a)The transition of a first excited state star to the ground state is shown. Here the radial metric is plotted against 
radius for various times starting from t = and then for ta = 250, tb = 500, tc = 1245, td = 4000, te = 5000, t/ = 5505, tg = 6000 
and th = 6370. The initial unperturbed and perturbed configurations arc shown. (The perturbed configuration has the lower 
second peak.) The initial mass of the star was M = .901mpi/m after perturbation. The radial metric initially has two peaks 
indicative of a 1-node configuration. As the star evolves and goes to the ground state one peak disappears. This can be seen 
in the curves from tc-th- Once in the ground state it oscillates and finally settles into a stable ground state configuration, (b) 
The transition of a 3 node configuration of mass 0.91m%i/m. This star loses enough mass during the course of its evolution to 
move to the ground state. 
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FIG. 9. (a) A comparison of the manner of black hole formation of four exctited state configurations. The first frame is an 
S branch star of central density cr{0) = 0.1 that tries to go to the ground state but fails to. As the central density increases 
decays to black holes occur on a shorter time scale. A plot of the collapse of a ?7 branch of central density (j(0) — 0.3 is shown 
in the next frame. Stars get more dispersive as one moves farther along the U branch. A star of central density 0.4 shown in 
the third frame has a decay time close to the previous one. Decay times then start to increase. A star of central density 0.5 
close to the M = Nm point (~ a(0) = 0.541) disperses to over ten times its radius before collapsing to a black hole, (b) The 
decay time to black holes is plotted as a function of central density, for one node configurations (first excited states), of boson 
stars without self coupling. Perturbations are only due to the finite differencing efi^ects of the numerical scheme. To make the 
comparisons meaningful the 95% mass radius (which is our definition for radius of a star) of every configuration considered was 
covered by the same number of grid points. Configurations for which the central density ct(0) < a2 move to the ground state. 
This value of the critical density a2 = 0.048 corresponds to a mass M = 0.91mp(/m. ai = 0.021 is the value of the central 
density corresponding to a mass M = 0.633mp; /m, which is the maximum mass of a ground state boson star. The decay time 
decreases with increasing cr(0) and this continues even for (7(0) > 0.25 which is the point of transition from the S-branch to the 
iZ-branch. The decay time then starts to increase as one approaches the M = N m point corresponding to a central density 
cr(0) = 0.541, beyond which the stars disperse to infinity rather than become black holes. The dispersion times of two such 
stars to ten times their original radius (95% mass radius) are also shown on the figure. 

FIG. 10. (a) The evolution of the metric function = —gtt for a A = 30, a(0) =0.1 boson star in the first excited state, with- 
out any explicit perturbation put, is shown. The configuration lies on the S-branch and has an initial mass M = 1.7A3m%i/m. 
The various time slices correspond to times t — 0,ta = 1060, tb = 1950 and tc = 1985. The lapse function collapses as an 
apparent horizon is approached, signaling the formation of the black hole (indicative of an inherent instability of excited states). 

(b) The evolution of the metric function = —gu for a A = 0, cr(0) = 0.122 boson star in the first excited state, without 
any explicit perturbation put, is shown. The configuration lies on the ^-branch and has an initial mass M = 1.23m%i/m. The 
various time slices correspond to times t = 0,ta = 450, it — 720, tc — 750 and td = 770. This star has a radius of 20.7 which is 
about the same as that of the configuration in (a). Again, the lapse function collapses when an apparent horizon is approached, 
as a black hole is being formed. The time scale of collapse is much less than for the A = 30 case in part (a). 

FIG. 11. (a) The initial field configuration of a 5- node star. The field has 5 nodes or extrema. The absolute value of each 
extremum is clearly smaller than the one preceding it. (b) A contour plot of a perturbed 5-node star, that ends in a black 
hole showing p x as a function of distance (vertical axis) and of time (horizontal axis) is shown. The density p is highest at 
the origin and has five other local maxima, each smaller than the previous one. The value of the maxima at the end are very 
small compared to the earlier ones, and to enhance the features p x rather than just p is shown in the plots. Each set of 
lines represents the maxima of p x and the number of lines in a set gives an indication of the height of the maximum. This 
particular configuration has a central density of a{0) = 0.075 and initial mass M = 3.07mp;/m. This star cascades through 
an intermediate 4 node state (around t = 1000) before proceeding to form a black hole. Cascades are characteristic of excited 
boson star decays similar to atoms in excited states going through intermediate states when transiting to the ground state. 

(c) The equilibrium density function of a 4 node star of central density a = 0.06 (right frame) is placed alongside the contour 
plot of the 5 node star described in 116 up to a time oi t = 1000 when it has gone into a 4 node state (left frame) . This plot 
shows that the transition of the 5 node star is to a perturbed 4 node state close to the one shown in 11c before it continues its 
evolution to a black hole. 



FIG. 12. The equilibrium profiles of a A = 800 star with central density a = 0.05 derived from the high A approximate 
equations and the exact ones are compared. A Schwarzschild exterior is attached to the approximate solution after the field 
vanishes. The three plots show the field a, Qrr and poo respectively. Clearly, the approximation matches the exact solution 
very well. 

FIG. 13. The mass of a high A star generated from the approximate equations is plotted as a function of a (cr/As). It 
shows the same basic structure as the profiles generated for low A using exact equations. The peak is at about .22A^/^mp;/m 
which means to achieve O.IM© would take A of order 10** a very large star to evolve numerically. Also plotted is N m {N is 
the particle number and m the mass of a boson). The crossing point of the two curves represents transition from negative to 
positive binding energy. 
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TABLE I. The ratio of the QNM frequency for A = 1600 to the QNM frequency for a given A is compared to (which is 
the predicted ratio for large A) for A = 1200, 800 and 1600. The higher A values match better as expected. The initial central 
density is ct(0) = 0.4. 
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